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An interface description and numerical simulations of model A kinetics are used for the first 
time to investigate the intra-surface kinetics of phase ordering on corrugated surfaces. Geo- 
metrical dynamical equations are derived for the domain interfaces. The dynamics is shown 
to depend strongly on the local Gaussian curvature of the surface, and can be fundamentally 
different from that in flat systems: dynamical scaling breaks down despite the persistence 
of the dominant interfacial undulation mode; growth laws are slower than t 1 ^ 2 and even 
logarithmic; a new very-late-stage regime appears characterized by extremely slow inter- 
face motion; finally, the zero-temperature fixed point no longer exists, leading to metastable 
states. Criteria for the existence of the latter are derived and discussed in the context of 
more complex systems. 
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1 Introduction 

Many two-dimensional surfaces exhibit internal degrees of freedom which allow for phase 

ordering or phase separation to occur within the static or dynamic surface. Examples are 

lipid bilayer membranes [1], crystal growth on curved surfaces [2], and thin film deposition [3]. 

Interaction between the shape and internal degrees of freedom of the surface are believed 
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to play an important role in such systems, initiating, modifying or eliminating chemical 
or physical intra-membrane domain-ordering processes. The dynamics can be induced not 
only by varying the temperature, but also by changing such variables as the pH or ionic 
concentration of an aqueous embedding solution (in the case of lipid bilayers), or the lattice 
mismatch and deposition rate (film growth). Experimentally [4], both spinodal decomposition 
in binary mixtures and phase ordering in one-component systems have been observed. Phase 
ordering is more relevant to non-fluid phases and, in the case of biological lipid membranes, 
is possibly involved in the control of enzyme adsorption and protein-enzyme interaction on 
the surface of a bilayer. 

For curved surfaces, the consequences of the interplay between surface shape and intra- 
surface pattern formation are at present largely unexplored in literature. The small amount 
of experimental results are rudimentary, and most of the theoretical work has been limited, 
due to the mathematical and numerical challenges involved, to equilibrium models and shape- 
perturbation calculations around equilibrium. Many of the theories assume that the domains 
are already formed (see, for instance, [5]). Typically, some form of bilinear coupling between 
an intra-membrane order parameter and the local geometry of the surface, such as mean 
curvature, is used. When kinetics are considered, one obtains a dynamical equation of 
state similar to that of the Random-Field Ising model with long-range correlations [6], but 
with the long-range correlations entering through the gradients of the intra-surface order- 
parameter [7]. One result has been the creation of shape phase-diagrams [1, 8]. Recently, some 
researchers have done simulations to investigate shape change in surfaces made of two types 
of lipids[9, 10]. These simulations, relying on Monte Carlo and bulk Langevin equations, 
are extremely computer intensive, limiting the possible complexity of the surfaces or the run 
time of the simulations to experimentally irrelevant cases. 

Any dynamics occurring in curved spaces invariably introduces new and non-trivial con- 
cepts and effects which require considerable care and exploration. The problem of how 
the dynamical shape of a surface and an intra-membrane ordering process may affect each 
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other is complicated. We therefore first consider the simplest case where there is no explicit 
coupling between the two, and new correlations arise through the geometry of the curved 
surface. Thus, rather than inquiring about the shape-change when the membrane is inhomo- 
geneous, we investigate the influence of non-euclidean geometry on pattern formation and 
phase ordering kinetics when it occurs on a static curved surface. By comparing the results 
with known results for model A kinetics within a flat (euclidean) surface, we identify the 
sole effects arising from surface curvature. We hope that this will subsequently allow for 
a systematic extension of the problem in which the surface is not static and the explicit 
coupling between the intra-membrane order parameter and the local geometry of the surface 
is also included. 

More explicitly, we report in this paper novel results of a ground-up study of relaxational 
pattern-formation occurring within static two-dimensional surfaces, using model A[ll] as 
a starting point, without any explicit coupling to the local geometry other than through 
diffusion[7]. In section 2, we setup the analytical framework for the intra-surface phase or- 
dering kinetics on curved surfaces. The bulk kinetics are described by the Non Euclidean 
Model A equation, Eq.(l). Systems with phase ordering instabilities, with homogeneous 
initial state and small random fluctuations, often evolve to form domains separated by sharp 
interfaces. The bulk part of the domains equilibrates rapidly and the interface width satu- 
rates within a very short time. Further time evolution involves only interface dynamics in 
which interfacial widths hardly change, but the total interfacial length decreases in order 
to minimise the system free energy. In order to take advantage of this type of kinetics, we 
deduce in section 2, a set of equations [Eqs. (2), (13) and (28)], which are equivalent to the 
bulk Non Euclidean Model A Eq.(l). In this interface description, local geodesic curvature 
of the interface, local velocity of the interface and total interfacial length L(t) play a central 
role. One of the novel quantities that we investigate in detail in this paper (apparently for 
the first time in literature) is the Geodesic Curvature Autocorrelation Function. Equation 
(2) is the generalization of the well-known Allen-Calm equation for intra-surface interface 
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dynamics in curved surface systems. Equations (13) and (28) couple to Eq.(2) due to the 
geometry of interfaces on a curved surface. 

In section 3, we explain how the numerical simulations were carried out: the interface 
description became pivotal to the feasibility of simulations, providing gains of 50 - 100 
in runtime with regards to the standard bulk description, for the system sizes that were 
explored. We also describe the checks on the numerical algorithm that we performed for 
special cases where analytic results are known. 

For all the results described in sections 4 to 7, we first integrated the bulk Non Euclidean 
Model A Eq.(l) for a short time to generate sharp interfaces on a variety of surfaces, extracted 
the interfaces, and then evolved them in time using three coupled interface equations. Three 
important quantities used in these sections are K, K g , and Kg, which must not be confused: 
the intra-surface interfaces can be characterized locally by their total interface curvature K or 

— * 

by their geodesic curvature K g , while a curved surface, on which those interfaces move, can 
be described locally by its Gauss curvature radius .Kg [12]. These simulation results discuss 
how Kq affects interface evolution and what are some of the new features of Non Euclidean 
Model A dynamics on corrugated surfaces, such as non-power-law domain growth, a new 
dynamical regime characterized by extremely slow interface motion, breakdown of dynamical 
scaling, time-dependence of various dynamical lengths not relevant in flat systems, and 
metastable interface configurations and activated hopping. Two measurements of interest for 
characterizing the domain morphology are introduced, namely the autocorrelation functions 
for K g and K. We give some quantitative and qualitative analytical explanations for these 
results while pointing to some others which require further study. 

2 Interface Equations 

It was shown in [7] that on curved surfaces, the equation for model-A kinetics (sometimes 
refered to also as the time-dependent Ginzburg-Landau equation, or TDGL for short) must 
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be written as 

M- 1 ^ = 0-0 3 + e 2 V L 2 B (1) 

where M and £ are positive phenomenological constants, t is time and is the intra-surface 
order parameter, while Vl B is the Laplace-Beltrami operator [13]. Some phenomenological 
parameters have been eliminated by appropriately rescaling space, time, and 0[14]. We 
refer to Eq. (1) as the Non-Euclidean Model A equation. The order parameter could be for 
instance the local magnetization at the surface of a corrugated anisotropic Ising ferromagnet, 
following a quench from a high-temperature, disordered (i.e. paramagnetic) state in thermal 
equilibrium, to a low-temperature, thermodynamically unstable region of the temperature- 
magnetization phase-diagram.. In this equation, there is no constraint on the average order 
parameter per unit area as a function of time. This differs from the well known model 
B which describes spinodal decomposition in binary mixtures [14], and for which the order 
parameter is conserved. From Eq. (1) we derived in [7] a Non-Euclidean Interface Velocity 
equation for interfaces on curved surfaces, 

V = Mi 2 K g (2) 

where z7 is the local interface velocity and K g is the local geodesic curvature of the inter- 
face on the surface. Equation (2) reduces to the well-known Allen-Calm equation in the 
Euclidean limit [15]. Note that the derivation of Eq. (2) given in [7], although concise, is not 
as geometrically transparent as the different one given in [16]. 

Further insight into the interface dynamics can be gained by considering the evolution of 
Kg in time. We give here a parameterization-invariant derivation of geometrical dynamical 
equations valid for any model-A interface on any Riemannian surface. To do this, it is 
necessary to parameterize the interface. 

Because, as indicated by Eq. (2), the interface moves normal to itself, it is convenient to 
use a parameterization of interfaces which exploits this feature: the normal gauge [17]. This 
dynamical constraint defines the gauge, whose units of length are therefore not physical 
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units, and change both in time and space. We denote the parameter of the normal gauge by 
a. Physical measurements, on the other hand, are always done using real (hence constant) 
units of distance. The parameterization which satisfies this is the arclength gauge, whose 
parameter we denote by s. 

The two gauges are linked by a metric which we denote ^fg. Denote the curvilinear 
coordinate system on the surface by (w^w 2 ) = (u,v). The surface is described by a series 
of three-dimensional vectors X = [x(u, v), y(u, v), z(u, v)], while the interface is described 
by a series of two-dimensional "vectors" R(a) = v(a)] or equivalently by a series of 

three-dimensional vectors R(a) = [x(a),y(a), z(a)]. We use the convention of an arrow to 
denote a three-dimensional vector, and a tilde a two-dimensional vector in the tangent space 
of the surface. A typical model-A interface, labelled J, is shown in Fig. 1, along with 
the vector R going from some arbitrary origin on the surface to a point P(a) on /, and the 
tangent f and normal n to the interface at P. Every point of an interface corresponds to a 
unique a which does not change in time, since the interface moves normal to itself. 

The relationship between the differential elements of length ds and da is given by 

ds = y/g da. (3) 

where the metric is defined by 



g(a,t) 



2 

dR 



(4) 



da 

Because da is constant in time while ds is not, g changes in time. Its evolution equation 
must be determined, making use of the requirement that a point of coordinate a moves only 
perpendicularly to the interface. This can be done in the following way. 
First note that f and V are defined as: 

f S 2£ = («>') (5) 
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where a prime denotes differentiation relative to the arclength s. Taking the time derivative 
at constant a on both sides of Eq. (4), we have 



dg\ 2 dR ( d_ OR 
dt /„ da \ dt\ da 



Now since 



, da 1 
« = 7TT = ^ ( 8 ) 



9s ^ 
d 2 q -1 
(9s 2 2g 2 da 



a = ttt = ^3^-' ( 9 ) 



the Laplacian in arclength s is 



d 2 1 d 2 d 



ds 2 g da 2 da 

so that from the definition of the (three-dimensional) interface curvature vector K = , 



ga2 9 {K -tcCrf). (11) 

Substituting Eq. (11) in Eq. (7), using the orthogonality between V and f and rewriting 
Eq. (7) in terms of ^fg rather than g one obtains 

« 12 » 

which further reduces to 

^f)=-^K t -V, (13) 

when one notes that V is parallel to n and lies within the surface, whereas K has a component 
parallel to V and one normal to the surface. I.e., if the normal to the surface is denoted by 
N, then K = \K g \n + \K n \N, with K n the three-dimensional curvature of a geodesic tangent 
to I at P. Equation (2) and the evolution equation for g, Eq. (13), are two of the three 
interface equations that are used in our interface description. 



Now we proceed with a similar but lengthier derivation for (dK g /dt) a . The dot product 
of n with Eq. (11) yields 

5 <9a 2 

Applying the time derivative and using the chain rule, 

d\K q \\ n d 2 ( dR\ n d 2 R(dg\ 1 f dh\ d 2 R , , 

1 1 ' ' :k ' a"?- (15) 



<9t / g da 2 \ dt I g 2 da 2 \ dt J g\dt J da 2 

/a \ / a 

The second term is simplified via Eq. (14) and Eq. (13). The third one is simplified by using 
Eq. (11) and noting that n can only change in direction, hence (ff ) is oriented along f. 
With those simplifications one gets 

d\K a \\ n d 2 V „, r -t, r -+ _ _ /<9n\ 

=7s? +2|A » |AV ^^ aV 'Ui (16) 

The second derivative of V is 

d 2 V d 2 , <9 2 n dhdv ^d 2 v , „. 

H = ^ + 2 7 r- 7 r- + n— -. (17) 



9a 2 9a 2 (9a 2 9a <9a 9a 

The first term is simplified by writing the second derivative of n in terms of arclength via 
Eq. (3), applying the chain rule and using the two Frenet equations for curves in Riemannian 
spaces [13] 

dt 

— = \K g \h = K g (18) 

| - -H*. W) 

This yields 

9^ = -^rg- s (\K g \^)-g\K g \K. (20) 

The second term of Eq. (17) vanishes when dot-multiplied with h given the second Frenet 
equation. After substitution of Eq. (17) and (20) in (16), we have 

d\K g \\ ld 2 v _ ^ ... (dn\ 
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Finally we show that y/gr ■ (||) in the last term of Eq. (21) satisfies 
It is easiest to show this by first deriving The norm of the velocity is 



v — n 



dR 
~dt 



Taking the arclength derivative and using the chain rule, 

dv n d I dR\ I dR\ dn 

ds = ^/f "da\dt) + \9t) "9s 

n ( d ( dR\\ 

since V is orthogonal to |j. Using ^ = §f and the chain rule again yields 

— (a 

where the second equality is obtained by using d(f -n)/dt = 0, completing the demonstration. 

Thus, substituting Eq. (23) into Eq. (21) and using Eq. (10), we get the last of the three 
curvature equations: 

/ a 

The three equations, Eq. (2), Eq. (13) and Eq. (28), are a coupled set of curvature equations 
describing the evolution of interfaces in the Non Euclidean Model A. Note that in these 
equations, if two-dimensional vectors are used, the scalar product must be properly defined 
via the surface's covariant metric tensor = J£ • |£, i.e. a ■ b = aig^bj = a % gijb> , where 
an implicit sum over alternate repeated indices is assumed. Equation (28) and Eq. (13) are 
purely geometric but only valid in the normal gauge and when z7 is parallel to n. 
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The first term on the right hand side of Eq. (28) is diffusive and causes all modes of K g 
to dissipate, except the zeroth mode (K g = 0) which is not affected by it. The diffusive 
term thus seeks to shrink interfaces either into straight lines (on a curved surface, geodesies) 
which have zero K g , or perfect (i.e. geodesic) circles which have constant K g . On the other 
hand, due to the Non-Euclidean Interface Velocity equation, the second term on the right 
hand side of Eq. (28) is cubic in K g . It seeks to increase the curvature and dominates 
when the diffusive term is negligible, i.e. for straight lines and circular domains. Hence it 
causes circular domains to shrink in radius, as is indeed the case in Euclidean model A. If 
the circular domain is also a geodesic, so that K g = everywhere along the interface, the 
domain will not change in shape. When the model-A dynamics proceeds from a quench of 
a system from a high temperature, such that the domains are initially very disordered, the 
diffusive term of Eq. (28) dominates, broadly speaking, during the early stages of interface 
motion, whereas the cubic term dominates only in interfaces which have become circular. 
However, this was verified in numerical simulations to be only approximately true. 

Another interesting aspect of the geometric equations is that substitution of Eq. (2) 
into Eq. (13) shows how ( is always negative. Therefore the distance between close- 



by interface points always decreases in time, and by extension also the total length of an 
interface, which is given by L = f y^da: 



Without these results it is difficult to ascertain whether Eq. (2) is consistent with the physical 
nature of the bulk equation for model A kinetics, where the quantity of order-parameter 
gradients (equivalent to the total length of interface in the system) should monotonically 
decrease in time. 

The apparent simplicity of the Non-Euclidean Interface Velocity equation can be mis- 
leading. The geodesic curvature of a line on a curved manifold depends both on the position 
of that line and its orientation on the surface. Therefore, K g introduces into the dynamics 
a new, position- and orientation-dependent length scale. This strongly suggests that Non 




a 




(29) 
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Euclidean Model A dynamics cannot be self-similar as it is in flat systems, and therefore 
that dynamical scaling will not be observed on curved surfaces, except perhaps on self-affinc 
surfaces, such as found over a certain range of length scales in lipid bilayer membranes. 
Different surfaces should show different phase-ordering kinetics. This is further evidenced 
in Eq. (29) where L — proportional to the reference length scale of the dynamics, if it is 
present [18] — depends on f \K g \ 2 ds. 

Nonetheless, a study of Non Euclidean Model A on the torus manifold, reported in [7], 
revealed no clear signature of the non-Euclidean nature of the surface when investigated 
through the evolution of the interface density (quantity of interface per unit area) and 
the A-structure factor [19] of the order parameter. The latter is a two-dimensional order 
parameter structure factor which depends only on the intrinsic geometry of the surface 
and the order-parameter configurations. This absence of signature may be due not only 
to the coarse-grained nature of the A-structure factor, but also to finite size effects, since 
domains feel the geometry only when their size becomes comparable to the torus size. This 
suggests that surfaces with geometrical features closer to those of model-A domains at early 
times should be investigated. Indeed, as we show in section 6, interfaces can get caught 
around certain surface bumps, causing a drastic slowing down and even immobilization of 
the model-A dynamics. Before investigating this, we outline the numerical methods used in 
the simulations. 

3 Numerical method 

It is customary to simulate equation for model A kinetics by discretizing them and evolving 
the order-parameter configurations via an Euler integration scheme. The main limitation of 
such discretization schemes is that the time step for evolving the system is limited by the 
smallest space mesh in the system. On curved surfaces, this can be — and most often is - 
a severe impediment, as the surface mesh can rarely be optimized so as to be homogeneous. 
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Doing this involves computational challenges in its own right. Slightly more sophisticated 
methods such as spectral, predictor-corrector and implicit schemes are possible but the 
stiffness of the model-A kinetic equation still induces strong dependencies on spatial mesh 
which make them unappealing. 

We have found that using Eq. (2) as the basis for an interface description has many 
advantages over a bulk description as given by a model A equation. Namely, the quantity 
of information to manipulate at every time step is an order of magnitude smaller at the 
beginning of the simulation than with a bulk description, and rapidly decreases in time, while 
for the bulk description it remains constant in time, even if the order-parameter configuration 
has become completely homogeneous. Secondly, the surface can be discretized independently 
from the interfaces so that numerical instabilities associated with the surface mesh disappear. 
Thirdly, the interfaces being one-dimensional, the integration algorithms remain fairly simple 
and more sophisticated algorithms are much easier to implement. Finally, certain quantities 
of interest, which we introduce below, are far easier and faster to compute from the interfaces 
than from the bulk configurations. 

One inevitable drawback of an interface description is that it cannot describe how in- 
terfaces form from the small order-parameter fluctuations present just after the quench. 
Therefore, the initial interface configurations must be obtained independently. The easiest 
way to do this is of course to use the bulk description to evolve the order parameter for only 
the relatively short time necessary for interfaces to become established, and then extract the 
interfaces from the resulting configurations. Even for short integration times, this suffers 
from the drawbacks mentioned above. We give one alternative, in appendix C of [16], devoid 
of any bulk steps and dependent only on the statistical properties of interfaces. 

For this reason, initial interface configurations are always obtained by integrating the 
bulk equation (1) from t = + to t = 17, when the domains have fully formed and interfaces 
can be extracted. For integration of bulk equations on curved surfaces, [20] has proved 
invaluable. Once the bulk configurations are obtained, the interface extraction is done with 
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a program which we refer to as geninterf. A full-length run, done by a program we refer to 
as IFinteg, then takes the interface configurations from t — 17 to t — 10000. Such long runs 
were necessary in Non Euclidean Model A due to the slowing down of the dynamics at later 
times. A typical initial interface configuration is shown in Fig. 2, as would be seen in flat 
systems and in surfaces whose geometric features are much larger than the size of the initial 
interface convolutions. 

An interface represented by its ordered set of coordinate points R(a) is evolved in the tan- 
gent space (u 1 , u 2 )(= (u, v)) of the surface via the Non-Euclidean Interface Velocity equation 
by writing the latter in the form 

fdu , , ... 2 



ym) = M ^ K ^ (so) 

where n l is the i th component of the interface normal n, and K g is given, in compact notation, 
by 

K g = eytiV" + n^'u 1 ') (31) 

where again implicit summation over repeated indices {i,j = 1,2} is assumed. The prime 
denotes differentiation with respect to arclength. Arclength is always the physical length of 
the curve as measured in the Euclidean embedding space of the surface, e^- is (j — i)^fg~^-, 
with g m = g\\gi2 — gngvi the determinant of the metric of the manifold rather than that of 
the interface, which is y/g. Equation (31) defines K g as a signed scalar, therefore n is defined 
as f rotated by 7r/2 counterclockwise and requiring n ■ f — 0, leading to 

n = —}=(-v'g22 - u'g 12 , u'gu + v'g 12 ). (32) 

The numerical integration of Eq. (30) for the j th interface point is done by first computing 
yfg~m~ and ^fg at the point, using these to obtain K g from Eq. (31) and n from Eq. (32), and 
then by moving the point via the time map 

u% t + At) = u l (j, t) + AtM?K g (j, t)n% t);i = 1, 2 (33) 
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This Euler scheme becomes unstable if the discrete time step At is too large. With the 
a parameter equal to the index of a point along the interface, ^fg is the physical distance 
between points. Therefore, each interface is evolved with its own time step At = 0.4At max , 
with 

At max = |(v^ min ) 2 (34) 

and y / ^ min is the smallest distance between two neighboring points of the discrete interface. 
As the integration proceeds, yfg min decreases, which forces a decrease of At. For this reason, 
the interface is remeshed with all points separated by a distance of 1 unit whenever y/g { < 
0.5. The remeshing algorithm can be made very efficient, one remeshing requiring barely 
more than a few integration steps. In flat systems, the interface description is roughly 5 
times faster than the bulk description. On curved surfaces, the gain increases dramatically: 
it was 50 or 100-fold in the simulations reported here, and will be much higher for more 
complex surfaces. Simulations using only the bulk description require, for a batch of 40 runs 
on the type of surfaces used here, on the order of 60 to 80 days on an HP9000s735. 

The integration algorithm was verified for correctness and accuracy through several in- 
dependent checks. Notably, the Non-Euclidean Interface Velocity equation can be solved 
exactly for a circular domain in a flat system. The simulation of such a domain yielded a 
curve indistinguishable from the theoretical prediction. Statistical runs starting from random 
initial order-parameter configurations were done with both bulk and interface descriptions 
and compared by measuring the amount of interface per unit area and comparing interface 
configurations and curvatures. Differences smaller than the interface width were found for 
the configurations, while the dynamical growth exponent of domains in flat systems was 
found to be 0.48 ± 0.01 with the interface description, closer to the theoretical prediction of 
1/2 than the bulk numerical integration result of 0.45 ± 0.02. The interface description is 
hence favored if only in accuracy. 

Simulations were also done on the torus manifold with both the bulk and interface de- 
scriptions and found to be once more identical apart from the similar difference in the 
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dynamical exponents. In general, the interface description provides more accurate results 
than the bulk description, as seen by repeating the simulation of a band domain on a torus 
manifold[7, 16] for different surface mesh coarseness. The interface velocity from the bulk 
description converges towards the Non-Euclidean Interface Velocity theoretical prediction 
as the surface mesh is refined, while that from the interface description falls exactly on the 
theoretical prediction for all surface meshes used. This is shown in Fig. 3 where v{t) is the 
interface velocity, and Qi the Interface position on the torus (the torus manifold with its 
coordinate system and manifold are shown below in Fig. 5). 

Finally, as a further accuracy check, a circular domain on a surface consisting of one 
circular bump was simulated as the interface equations can be solved analytically. With the 
surface defined by 



with A the amplitude of the bump and a its half-width, the relationship between K g and 
the radius R of the domain is 



The comparison of the theoretical prediction and numerical simulation is shown in Fig. 4 
where the product K g R is plotted vs R. Two different surface meshes were used, showing 
what kind of surface coarseness is sufficient for accurate results: the longest surface tether 
should be around 0.7 units long for best results, though 1.3 units gave satisfactory results as 
well, as seen in the figure. 

4 Effects of surface curvature on domain growth 

A naive interpretation of the Non-Euclidean Interface Velocity equation may suggest that 
Non Euclidean Model A dynamics should always be slower than Euclidean dynamics: \K g \ is 

— * 

always less than or equal to the total interface curvature \K\. However diffusion occurs faster 



X= U) v,Ae- {u2+v2)/2 > 



(35) 



1 



(36) 




15 



(slower) in regions of negative (positive) surface Gaussian curvature Kg [12], because there 
is more (less) area available for a given interface length in a region of negative (positive) 
i^G [21]. This suggests that interfaces arising from the Non Euclidean Model A equation 
should disappear more slowly where Kq > 0, but faster where K G < 0. Yet, simulations 
of the Non Euclidean Model A on the torus manifold, starting from random initial order- 
parameter configurations, showed no dependence whatsoever on Gaussian curvature [7]. This 
may have been due to the use of the A-structure factor [19] as the measure of the two- 
dimensional Order Parameter Structure Factor or to the combination of the topology of the 
torus with the percolating domains. 

To settle the matter, one ovoid interface was therefore simulated in the K G > region of 
the torus manifold, another in the Kq < region and again in a flat system (Kq = 0). The 
torus manifold is convenient for many reasons discussed in [7], most importantly that two 
well-separated regions of oppositely signed Gaussian curvature can be defined on the torus. 
A drawing of part of a torus surface with coordinates and parameters is shown in Fig. 5. 
The small radius of the torus is denoted r and the large one R. The initial radius of the 
ovoid is such that it is completely included in the appropriate region of the torus, i.e. nr /2 
as measured in (6>, ip) space (see Fig. 5). 

Figure 6 shows a log-log plot of — dL/dt as a function of l/L, where L(t) is the time- 
dependent length of the ovoid interface in the system. A straight line is obtained in the 
flat system, as expected for Euclidean model A. However, differences as large as a factor of 
two can be seen at early times with the curved regions, when the domains are large (roughly 
20% of system area). The difference remains substantial until the domains are very small. 
This suggests that the absence of visible non-Euclidean effects reported in [7] could be a 
characteristic of the A-structure factor. Note that the oscillations for the K G ^ curves are 
due to the numerical discretization of the torus manifold. 

The above simulation result can be understood analytically with the help of Eq. (13) 
by considering surfaces of constant Gaussian curvature (i.e., a sphere, a flat plane, and a 
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hyperbolic plane). In this case, a circular domain remains circular, as indicated by Eq. (28), 
i.e. the geodesic curvature is constant everywhere along the interface at any given time 
(but changes in time). Integrating both sides of Eq. (13) along the interface then yields 



in a curved surface should be of the form ^Kq, with 7 a dimensionless geometrical factor 
of the order of 1. Moreover, when Kq is increased from (i.e. the surface is a sphere of 
decreasing radius) while keeping L constant, the geodesic curvature monotonically decreases, 
until it eventually becomes when K G = (2-k/L) 2 . On the other hand, when K G is increased 
negatively from (i.e. the surface is a hyperbolic plane), the area available to the circular 
domain increases, such that in order to keep L constant, the radius of the domain must be 
decreased, increasing the geodesic curvature. This leads to the first order approximation for 
the time change of the circular domain's interface length, 



with 7 > 0, providing a qualitative explanation the difference in growth laws for the ovoid 
domains in regions of different Gaussian curvature. 

Given the physical origin of the influence of the surface Gauss curvature Kq on domain 
growth in the Non Euclidean Model A, much the same can be expected in non-Euclidean 
model B as well as in other systems where diffusion and interfaces are present on curved 
surfaces. This slowing down for positive Kq was not observed in numerical simulations of 
pure model B on the static sphere [9], because of the short run times used. However, similar 
(though not identical) Gauss curvature effects were seen in simulations of crystal growth on 
toroidal geometries [2]. In the context of phase-ordering kinetics and phase-separation in lipid 
bilayers, where diffusion is known to play a very important role biologically, this dependence 
on Kq could be of use in controlling the rate of such processes as protein diffusion via 
membrane shape change, or inducing enzyme-protein interaction in limited parts of a cell. 



— dL/dt ~ LKg, where K g is a function of time (note also Eq. (29)). It is useful to express 
Kg in terms of Kq. In a flat system K 2 = An 2 / L 2 when the domain boundary has length 
L. For small Kq 7^ 0, dimensional analysis indicates the first order correction to K 2 when 




(37) 
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5 Simulations on sinusoid surfaces 



When the initial order-parameter configuration is one of complete disorder, the non-linear 
regime of Euclidean model A is characterized by the relatively slow motion of sharp, con- 
voluted interfaces delimiting domains of = ±1 phases. There are two characteristics 
of Euclidean-model-A dynamics which are particularly important here. First is the self- 
similarity of the dynamics, leading to dynamical scaling: system configurations at a time t\ 
look statistically identical to configurations at an earlier time to, if they are rescaled length- 
wise by an appropriate factor. Hence all dynamical lengths in the system have the same time 
dependence, so that they can all be expressed in terms of one arbitrarily chosen reference 
length scale L(t). The numerical value of L is not as important as its time dependence, 
which is the second characteristic of relevance: L(t) ~ t 1 / 2 . 

It is common to refer to L as a dominant or typical length scale in the dynamics, but 
the order parameter structure factor for Euclidean-model-A systems shows a peak at zero 
wavenumber, indicating model-A dynamics does not have a dominant length scale, only a 
unique time dependenceior all dynamical lengths. We have shown for the first time [16, 18], by 
considering the curvature correlations along the interface, that Euclidean-model-A systems 
exhibit a dominant undulation mode not present in the order-parameter structure factor. 
If one is to talk of a dominant dynamical length in model A, it is the wavelength of this 
undulation, not the width of the order parameter structure factor. This suggests it has 
a different nature from that of Euclidean model B, where the dynamical length is visible 
in the order parameter structure factor. Notably, the absence of the dominant dynamical 
length in model A's order parameter structure factor may be due to the absence of any phase 
information in it. 

Now consider a surface consisting of a large array of bumps: 



Experimentally, fluid as well as crystalline lipid bilayers are known to adopt a similar shape 
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(called "egg carton") under certain conditions [22]. Also, more complex surfaces such as self- 
afflne surfaces common to fluctuating lipid bilayer membranes share many of the qualitative 
features of sinusoid surfaces. 

Simulations of Non Euclidean Model A on sinusoid surfaces, starting with random initial 
order-parameter configurations, were done for several values of A and A, but here we focus 
on the A = 20, A/\ = 0.2 surface, with system sizes of 100 x 100 to 300 x 300. 

The normalized geodesic curvature autocorrelation function is defined by 



where (...) denotes an ensemble average with an average over all interface points s . It is 
plotted in Fig. 7, where the horizontal axis has been rescaled with L(t), defined here as 
first zero of the geodesic curvature autocorrelation function, and the error bars (not shown 
for clarity [23]) are much smaller than the vertical offsets of the curves. In flat systems, as 
explained in [16, 18], all such geodesic curvature autocorrelation functions fall on top of one 
another, due to dynamical scaling. Figure 7 thus shows that for the non-Euclidean case, 
dynamical scaling breaks down, as expected from our discussion of Eq. (2) but contrary to 
runs on the torus manifold. The other important feature is that as time increases, the domi- 
nant interface undulation mode monotonically decreases in intensity, signifying the system's 
degree of order decreases as the interfaces explore ever larger length scales. This is sensible 
given that when the dominant length is much smaller than the geometrical features of the 
surface, the interface correlations decay fast enough to become negligible. As the undulation 
wavelength increases, the correlations become negligible at increasing distances, introducing 
geometrical variations in the interface shape, thus increasing the degree of disorder in the 
system. 

These two features — breakdown of dynamical scaling and existence of a dominant dy- 
namical length — further confirm that the presence of a dominant dynamical length does 
not guaranty dynamical scaling. 




(K g (s ,t)K g (s + s,t)) 



(39) 
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Figure 7 also indicates that the first zero of the geodesic curvature autocorrelation func- 
tion can no longer be used as a reference dynamical length. Two well-defined dynamical 
lengths valid for interfaces on curved surfaces are the width a of the Gaussian part of 
the geodesic curvature autocorrelation function, and the inverse interface density. Though 
G Kg (s, t) is found to remain Gaussian on short length scales, the time dependence of its width 
is no longer a t x l 2 power law, contrary to flat systems. For the sinusoid surface used here, 
it was found to be logarithmic. The logarithmic law is shown in the inset graph to Fig. 7, 
where linear regression gives a 1 ~ (63 ±3) log 10 (t/to), with to = 6 ±2. This logarithmic time 
dependence is not a universal feature, as sinusoid surfaces with smaller values of A/ A showed 
power-law growth but slower than t 1 ! 2 . It seems likely that a sub- logarithmic growth law will 
be seen for larger A/X, similarly to results reported in [24] on the related Random-Field Ising 
Model. The width a can be used as the reference length scale L(t) to rescale the geodesic 
curvature autocorrelation function, as shown in Fig. 8. The inset shows the result for flat 
systems, exhibiting dynamical scaling, within error bars. The main graph shows the result 
for the sinusoid surface. There the curves superpose only on short length scales, while at 
long length scales they systematically become wider, with the drift being much larger than 
the error bars. 

While the width of the geodesic curvature autocorrelation function is a short-length-scale 
characteristic of the dynamics, the interface density is more sensitive to large-scale features. 
It is defined as l(t) = L(t)/A s , where L(t) is the total length of interfaces at time t and A s is 
the system area (constant in time). Measurements of l(t) for the sinusoid surface were done 
and compared with those in flat systems of various sizes. This is shown in Fig. 9, where 
—d(lnl)/dt is plotted as a function of l// 2 . This scale produces one unique curve for all flat 
systems, independent of system size or quantity of interface. Any straight line on this plot 
indicates a power law in time for £(t). If the slope of the line is denoted m, the intercept b 
and £o = £(0), then straightforward integration yields 

£(t) = (£ 2m -2ml0Y)^, (40) 
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indicating I ~ t 1 / 2 ™ a t late times. More negative m thus corresponds to slower interface 
dynamics. The interface density has units of one over length, while the dominant length 
scale in Euclidean model A grows as t 1 ^ 2 . Hence the scaling assumption valid for flat systems 
predicts £ ~ t~ 1//2 , i.e. m should be equal to -1 in flat systems. The linear regressions for 
the flat system curve (where bulk integration was used rather than the interface description) 
gives m = —1.07 ± 0.01, which means a power law of —0.46 ± 0.01 instead of -1/2. 

At early times (small values of l/£ 2 ), model-A dynamics on the sinusoid surface also 
follows the same power law. At later times however, one can distinguish a fairly long time 
regime, extending from t ~ 250 to t ~ 1600, almost a full decade of time, during which m 
is more negative. The dynamics has therefore substantially slowed down. Linear regression 
in this time domain gives m = —1.70 ± 0.04, corresponding to a power law behavior for I of 
t -0 - 3 . A batch of runs for a smaller system of size 100 x 100 gave m = —2.0 ±0.1, suggesting 
the statistical error on the slope could be as large as 20 to 30%. 

A very interesting characteristic of this plot is the presence of a very late time regime, 
where the dynamics is extremely slow, with the slope m = — 8 ± 2. This starts at t ~ 2000, 
independent of system size and discretization. Real-time animation of moving interfaces in 
this regime shows that the geodesic curvature is zero almost everywhere, but not in suffi- 
ciently many places to completely halt the dynamics. The interfaces, which waver between 
the bumps, are very slowly hopping the bumps one by one. 

6 Metastable states 

Corrugated surfaces (such as given by Eq. (38)) allow for local minima to appear in the 
configuration space of interfaces and therefore to trap the latter in metastable states, because 
Eq. (29) implies that Non Euclidean Model A interfaces can only decrease their length. 
The simplest example of this would be an ovoid interface circumventing two bumps on the 
sinusoid surface, as pictured in Fig. 10. Consider the interface as it shrinks in length. It is 
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forced to tilt itself, slowly moving on the outside towards the extremum of each bump. The 
total curvature at a point of the interface is the curvature as measured in three-dimensional 
embedding space. Recall that the total curvature vector if of a geodesic line on the surface 
is normal to the surface since K g = for a geodesic. Therefore a necessary condition for the 
interface to become stationary is for a configuration to exist where K, near both ends of the 
ovoid domain, is normal to the surface. A side view (x — z projection) of the bumps with 
the interface is pictured in Fig. 11 and a top view (x — y projection) in Fig. 12. This is not 
a sufficient condition however since all points along the interface must satisfy this criterion, 
but it does provide a minimal condition. 

We now denote the vector going from O to P in Fig. 11 as R — [x,Asm(kx)]. The 
tangent to the surface in P is therefore f = [1, Akcos(kx)}, and the condition is R ■ r — 0, 
leading to the transcendental equation 2x = —A 2 ksm(2kx). With the substitution x' = 2kx 
and 7 = Ak this is written 

Y x' 

This equation has a solution only when 7 is larger than a certain value, for when 7^0 
the left hand side tends to —00 while the right hand side is bounded between 1 and roughly 
—0.22. The condition is therefore — 1/7 2 > —0.22, or 

y > 0.34. (42) 

This was tested numerically for fixed A = 40. Simulations give a threshold of 0.42 ±0.02. 
This is an error of 20%, surprisingly good considering that the approximation uses a one- 
dimensional projection of the interface on the surface: it is not unusual for the presence of 
a second dimension to give a system more freedom, thereby softening constraints such as 
Eq. (42). 

It is important to note that Eq. (42) is a necessary (minimal) condition for interfaces to 
get blocked around two bumps, and that blockage around a larger number of aligned bumps 
occurs at similar or larger values of A/ A. This can be seen by considering three aligned 
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bumps instead of two. The middle bump constrains the interface to belong to the xy plane 
in its vicinity, so that point O of Fig. 11 still lies between bumps rather than, say, at the 
center of the middle bump. The condition Eq. (42) therefore still holds for any number of 
aligned bumps larger than or equal to two. This is sensible since what matters is A, not the 
number of bumps. 

However the situation is the reverse for bumps that are not aligned, which is more relevant 
with regards to the results of the Non Euclidean Model A simulations on sinusoid surfaces. 
Consider 4 bumps in a 2 x 2 array, and an interface that circumnavigates the four bumps. 
A top projectional view is shown in Fig. 13. Point O, one of the interface pivot points on 
segment AB, is still halfway between bumps, near the system boundary. The maximum- 
curvature director however no longer lies along x but along the line x = y. Therefore, the 
pivot point normal to the interface at P is at O' rather than O, decreasing the apparent 
wavelength by a factor of y/2, so that smaller amplitudes (by the same factor) are sufficient 
to trap the interface. In this case the threshold decreases to 0.24, with simulations giving 
0.270±0.006. The proximity to this threshold of A/A = 0.2 in the sinusoid simulations of the 
last section explains the existence of the very-late-stage, extremely-slow regime evidenced 
on Fig. 9. 

The same kind of argument can be extended to, say, 6 bumps in a hexagonal configuration, 
bringing O' yet closer to P. For a large surface consisting of a great number of bumps, 
gradually increasing A from zero should cause interfaces to become stationary first around 
very large conglomerates of bumps. However the larger the conglomerate must be, the 
rarer the occurrence. As A is increased further, smaller conglomerates of bumps trap closed 
interfaces, until an A is reached where the domains stop growing when the dominant interface 
undulation length becomes comparable to the surface A, leading to long-range disorder. 

The metastability thresholds can be generalized to more complex surfaces. For instance, 
lipid bilayer membranes are usually characterized by their bending rigidity k. Such mem- 
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branes form random self-affine surfaces whose average square width W 2 is given by[25] 

W 2 = ^f- (43) 

where L s is the size of the membrane as projected on the xy plane. The amplitude A can be 
approximated by the width W for a given size, while the wavelength A can be approximated 
by L s . As a consequence of Eq. (42), phase ordering may not occur at all if k < /cbT/8, as in 
this case Eq. (42) is satisfied on all length scales L s . As k is decreased towards /c#T/8, the 
Non Euclidean Model A dynamics should gradually slow down. When k < ksT/8, domains 
might still form but freeze in their early disordered state. 

For real systems the hopping condition on k is not likely to be as simple, since for such 
small values of k the surface tension a will usually be non negligible. The main difference 
is that the threshold involves k, a and the length scale of interest, without introducing any 
new difficulty This suggests that domains could order for some time, until their dominant 
interface undulation length becomes larger than a certain threshold value. 

In flat systems, it was shown by Bray [26] that model A exhibits a zero-temperature 
strong-coupling fixed point. The existence of metastable interface configurations on cor- 
rugated surfaces implies that thermal noise changes the quality of the dynamics once the 
domain interface undulations are on the same scale as the surface corrugations, thereby 
eliminating the fixed point. Consequently, in the very-late-stage (i.e. extremely slow) regime 
observed on sinusoid surfaces, interfaces will move predominantly via thermally-activated 
hopping. Similar conclusions should be valid for model B. Effect of thermal noise needs to 
be carefully considered. 



7 Total curvature correlations 

There is a second form of curvature autocorrelation involving the total curvature K of the 
interface, i.e. the curvature of the interface measured in three-dimensional space rather than 
in the two-dimensional space of the manifold. In three dimensions, a curvature scalar can not 
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be well defined. The vector of curvature is the only way to properly represent the curvature 
of the interface, so a dot product must be used in this case: 

&(.,«) 5 <*<°-J> •*<'■'»■ (44) 

(K(o,ty-) 

— * 

Note that for flat systems, where the scalar product of two K g is straightforward, both 
definitions of curvature autocorrelation, the total and the geodesic one, give the same func- 
tion. It is only in curved surfaces that they differ in very important ways. In the curved case, 
when an interface is stationary, it is so only because K g is everywhere and the interface is 
perfectly autocorrelated, from a two-dimensional point of view. This perfect autocorrelation 
stems from the fact that if two-dimensional observers within the manifold knew the position 
and orientation of the interface at one point of the interface, they would know the complete 
shape of the interface by solving the equations defining geodesics[13] with the appropriate 
initial values. The only disorder remaining in the system is that arising from the relative 
position of the domains and from the extrinsic geometry of the geodesic lines on the surface. 
This is not measurable via a geodesic curvature autocorrelation function. 

On the other hand, the total curvature of a convoluted interface at early times, when 
the domains are much smaller than any geometric length scale of the surface, is equal to 
the geodesic curvature of the interface, in modulus. Under this condition the total and 
geodesic curvature autocorrelation function are nearly identical at early times (analytically 
they are rigorously equal). At late times, when an interface has become stationary, the total 
curvature at a point of the interface is the curvature of the surface along the direction of 
that interface. This total curvature autocorrelation function will therefore give information 
about the curvature autocorrelations of the surface along geodesic lines. This information 
is part of the extrinsic geometry of the surface and does not enter the geodesic curvature 
autocorrelation function. The two autocorrelation functions are therefore complementary 
for interface dynamics on curved surfaces. 

When the interface correlations decay faster in space than the wavelength of the surface, 
Gg(s,t) is, as discussed earlier, qualitatively very close to G K (s,t). Both correlation func- 
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tions are not exactly equal as the scalar product in Gg(s,t) produces a slightly stronger dip 
amplitude by about 20%, and the error bars are substantially smaller in Gg(s,t). Near the 
minimum they are as much as 3 times smaller. However, while the average K g decreases in 

— * 

time, the K eventually starts increasing until it reaches values compatible with the surface 
A. Therefore, Gg(s,t) becomes very different from G Kg (s,t) at late times. 

Figure 14 shows G^(s,t) for the runs on the sinusoid surface, Eq. (38) with A = 4 and 
A = 20, at times between 17 and 10000 (17, 35, 50, 100, 150... 400, 500, 700, 1000, 1200, 
1450, 1750, 2100, 2500, 3000, 4000... 10000). The curves for t = 17 and 10000 are labelled, 
with curves at intermediate times moving gradually and monotonically from the former to 
the latter. For this surface, the wavelength of the bumps is not much larger than the initial 
dominant length scale of the interfaces in Euclidean systems, so that a cross-over regime is 
not seen. The correlations simply increase in time, as more and more interfaces coarsen while 
being pushed between the bumps and adopting the surface's wavelength. At late times, the 
interfaces can be approximated by a random sequence of straight lines of length A/2 and arcs 
of radius A/4 and arclength (7r/4)A/2 . This explains the peak positions in G^(s,t) being 
roughly at integer multiples of (tc/A)\/2. At early times, only two peaks can be seen, with 
a third peak starting to appear. The emergence of the peaks at larger distances seems to 
coincide with the Euclidean regime in Fig. 9. During the slow regime (250 < t < 1600), no 
more peaks appear. Those that have formed grow slowly in amplitude and saturate. Figure 
15 shows Gpc(s,t) for a sinusoid surface consisting of two wavelengths rather than one: 

X = [x, y, 4 sin(27rA;i;r) sin(27rA;i?/) + 6.25 sin(27r/c 2 a;) sin(27r/c 2 y)] (45) 

where Ai = 2n/ki = 40 and A 2 = 2ir/k 2 = 50. Hence both have A/X ~ 0.1 and the 
undulations of this surface are of longer wavelength than the previous one. In this case, 
the earliest curve looks exactly like G^(s,t) in flat systems and even has the dip positioned 
at roughly the same value, namely s ~ 10. Already at t = 24 the interfaces are being pulled 
into the valleys between the bumps. This is the crossover regime. The curvature correlations 
settle into their definitive value around t = 100, close to (tt/A)X/2 with A around 40. Much 
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longer wavelengths would be necessary to distinctly show the three regimes: a scaling regime 
at early times during which the amplitude of the minimum would not change, then the cross- 
over during which the dominant length scale saturates to a dominant mode of the surface, 
and finally the late-time regime where the interfaces move more slowly and hop the bumps. 
Note that contrary to the A = 4, A = 20 case, the very-late-stage regime of Fig. 9 never 
occurs because A/\ is not large enough. 

8 Summary and conclusions 

In summary, the non-Euclidean model A was used as a starting point for the more com- 
plex, Random-Field Ising Model like [6, 27, 24] models used in some lipid membrane [28] and 
related surface problems. By deriving a set of geometric dynamical equations, using the non- 
Euclidean interface velocity equation V = M^ 2 K g and developing an interface description, 
numerical simulations of this difficult system could be carried out and the results analyzed 
and understood qualitatively. Moreover, the techniques should be applicable without major 
impediments to dynamical surfaces and model B on curved surfaces. 

Whereas Random-Field Ising Model systems exhibit decelerated growth in the form of 
logarithmic growth laws of various kinds, here we have shown how model A on curved surfaces 
exhibit similar richness of dynamics without a bilinear coupling to the surface: (i) the time 
dependence of the amount of interface is strongly affected by the local Gaussian curvature 
of the surface, (ii) different dynamical lengths have different and even logarithmic time- 
dependencies leading to the breakdown of dynamical scaling without the disappearance of 
the dominant interface undulation mode, and (iii) metastable states exist above a threshold 
value of Aj A, leading to thermally activated hopping in the new very-late-stage regime. 
A more systematic study of the effect of K G , a better quantitative understanding of the 
observed time dependencies of L, and experimental methods to obtain the dominant interface 
undulation mode, would be useful at this point. 
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Figure 1: An interface / in the (u,v) space of the surface. A point P of the interface has 
parameter a and is given by a vector R(P) (bold R in figure), f and h are the tangent and 
normal vectors to I at P. 
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Figure 2: A typical early-time configuration of interfaces in a flat system. 
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Figure 3: Interface velocity of stripe domain as a function of Qj on the torus; R = 40, r = 20. 
Numerical simulation with A9 = 0.5 (circles) and 1 (squares), and analytical prediction (line) 
are shown. At t — both top and bottom interfaces are near the outside equator of the 
torus. 
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Figure 4: K g R of a circular domain as a function of its radius R, as given by numerical 
integration on 2 different grid sizes; also shown is analytical prediction eq. (36). Error bars 
too small to be seen. 
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Figure 5: Torus manifold of long radius R and small radius r$. Point P on manifold has 
coordinates (m 1 ,^ 2 ) = (0,ip). 9 = and 1 are at inner equator (labelled "ie"), 6 = 7rr is at 
the outer equator (labelled "oe"). 
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Figure 6: Log-log plot of — dL/dt vs l/L. Straight line is for simulation of flat system. The 
scatter for nonzero Kq is due to numerical discretization of the surface. 
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Figure 7: G Kg (s,t) at various times from t = 35 (long dashed) to t — 400 (short dashed), vs 
s/L(t). Sinusoid surface used is eq. (38) with A = 4, A = 20. Inset: Variance of Gaussian 
part of G Kg (s,t) vs time (times between 17 and 700). 
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Figure 8: Main graph: Same as main graph of Fig. 7 but using L(t) = a(t) to rescale 
G K (s,t). Inset: Same as main graph but for flat systems. 
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Figure 9: Logarithmic plot of dln(/)/dt as a function of l/l 2 , for sinusoid surface A = 4 and 
A = 20 (circles). Points for flat system are triangles and squares, lines are linear regressions. 
Error bars (not shown) are approximately same size as symbols. 



40 




Figure 10: Three-dimensional view of the 2-bumps surface with the ovoid interface (white 
line). System size 40 x 20, A = 10, A = 40. Surface is shown as transparent to interface. 
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Figure 12: Top projectional view for Fig. 11. Interface is thick line. K is curvature vector, 
other symbols correspond to those in Fig. 11. Line AB is "pivot" line for interface. Circles 
are contours for the surface. 
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Figure 15: Gg(s,t) for the bimodal sinusoid surface, eq. (45) with Ai = 40 and A 2 = 
Times are t = 17, 24, 35, 60, 100 with the remaining times at same values as in Fig. 14. 
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